library(sandwich)
library(mgcv)
library(tikzDevice)
library(reshape)
library(plyr)
library(survey)
source("panel-utils.R")
load("sensitivity-data-dem.RData")

#### Sensitivity analysis ####
sims <- 100
confound <- seq(-4, 4, by = 0.2)
sens.sha <- matrix(NA, nrow = length(confound), ncol =  sims)
sens.inc <- matrix(NA, nrow = length(confound), ncol = sims)
sens.win <- matrix(NA, nrow = length(confound), ncol = sims)
sens.winc <- matrix(NA, nrow = length(confound), ncol = sims)
sens.early.sha <- sens.sha
sens.early.inc <- sens.inc


strategic.bias <- function(alpha, a) return(alpha)
confirmation.bias <- function(alpha, a) {alpha * (2 * a - 1)}

dem.panel$y.a <- NULL
for (i in 1:length(confound)) {
  cat(i, "\n")
  this.alpha <- confound[i]
  dem.panel <- ddply(dem.panel, "race", transform,
                     tmp = sens.adj(y=demprcnt, a=d.gone.neg,
                       a.fits=dem.pscore, alpha=this.alpha, confound.func = confirmation.bias))
  newNames <- names(dem.panel)
  names(newNames) <- names(dem.panel)
  newNames["tmp"] <- paste("y.a.", i, sep = "")
  names(dem.panel) <- newNames
}


dem.final <- dem.panel[dem.panel$week == -1, ]
dem.final <- dem.final[!is.na(dem.final$sw),]




race.list <- unique(dem.final$race[!is.na(dem.final$sw)])
race.lengths <- c()
race.rows <- list()
for (i in 1: length(race.list)) {
  race.rows[[i]] <- rownames(dem.panel[which(dem.panel$race
                                             ==
                                             race.list[i]),])
  race.lengths <- c(race.lengths, nrow(dem.panel[which(dem.panel$race ==
                                                       race.list[i]),]))

}
names(race.rows) <- race.list
names(race.lengths) <- race.list

for (i in 1:sims) {
  cat(i, " ")
  if (i %% 10 == 0) cat("\n")
  races <- sample(race.list, replace = TRUE,
                  size = length(race.list))
  races <- as.character(races)
  star <- unlist(race.rows[races])
  bootid <- rep(1:length(races), times = race.lengths[races])
  dem.star <- dem.panel[star,]
  dem.star$bootid <- bootid

  dem.wmod.ni <- iptw(denominator = n.denom, numerator = numer,
                      data = dem.star, id = "bootid", time = "week",
                      family = "binomial", subset = deminc == 0,
                      ones  = !dem.star$pos.prob.subset,
                      bal.covs = strsplit(covs, "\\+")[[1]])

  dem.wmod.inc <- iptw(denominator = i.denom, numerator = numer,
                       data = dem.star, id = "bootid", time = "week",
                       family = "binomial", subset = deminc == 1,
                       ones  = !dem.star$pos.prob.subset,
                       bal.covs = strsplit(covs, "\\+")[[1]])
  dem.star$sw[dem.star$deminc == 1] <- dem.wmod.inc$sw
  dem.star$sw[dem.star$deminc == 0] <- dem.wmod.ni$sw
  dem.star$dem.pscore[dem.star$deminc == 1] <- dem.wmod.inc$d.pscore
  dem.star$dem.pscore[dem.star$deminc == 0] <- dem.wmod.ni$d.pscore

  dem.star <- dem.star[dem.star$week == -1,]
  dem.star <- dem.star[!is.na(dem.star$sw),]

  causal.des.star <- svydesign(ids = ~ race, weights = ~ sw,
                               data = dem.star)


  for (ii in 1:length(confound)) {

    share.form <- paste("y.a.", ii, " ~ ",
                        "deminc * d.neg.rec + deminc * I(d.neg.dur - d.neg.rec) + ",
                        base.vars, sep = "")

    share.form <- as.formula(share.form)


    inc.form <- paste("y.a.", ii, " ~ ",
                      "I(deminc == 0) * d.neg.rec + I(deminc == 0) * I(d.neg.dur - d.neg.rec) +  ",
                      base.vars, sep = "")

    inc.form <- as.formula(inc.form)


    causal.sha.star <- svyglm(share.form, design = causal.des.star)
    causal.inc.star <- svyglm(inc.form, design = causal.des.star)



    sens.sha[ii, i] <- causal.sha.star$coef["d.neg.rec"]
    sens.inc[ii, i] <- causal.inc.star$coef["d.neg.rec"]
    sens.early.sha[ii,i] <- causal.sha.star$coef["I(d.neg.dur - d.neg.rec)"]
    sens.early.inc[ii,i] <-
    causal.inc.star$coef["I(d.neg.dur - d.neg.rec)"]
  }


}

save(list = c("sens.sha", "sens.inc", "confound"),
     file = "sensitivity-output-dem.RData")

noninc <- data.frame(cbind(coef=sens.sha, se =NA, CI=sens.sha.ci,
                           alpha = confound))
inc <- data.frame(cbind(coef=sens.inc, se =NA, CI=sens.inc.ci, alpha =
                           confound))

load("sensitivity-output-dem.RData")

regTable <- each(mean, sd, lower.95.bound, upper.95.bound)
noninc <- adply(sens.sha, 1, regTable)[,-1]
inc <- adply(sens.inc, 1, regTable)[,-1]
colnames(noninc) <- c("Mean", "SE", "95% Lower Bound",
                      "95% Upper Bound")
colnames(inc) <- c("Mean", "SE", "95% Lower Bound",
                   "95% Upper Bound")

noninc$pos <- noninc[,3] > 0 & noninc[,4] > 0
noninc$non <- noninc[,3] < 0 & noninc[,4] > 0
noninc$neg <- noninc[,3] < 0 & noninc[,4] < 0
noninc.pos <- c(which(noninc$pos), max(which(noninc$pos)) + 1)
noninc.neg <- c(min(which(noninc$neg))-1, which(noninc$neg))
inc$pos <- inc[,3] > 0 & inc[,4] > 0
inc$non <- inc[,3] < 0 & inc[,4] > 0
inc$neg <- inc[,3] < 0 & inc[,4] < 0
inc.pos <- c(which(inc$pos), max(which(inc$pos)) + 1)
inc.neg <- c(min(which(inc$neg))-1, which(inc$neg))

### Sensitivity Figure ###

##############
## Figure 7 ##
##############

par(mfrow = c(1,2), cex.main = .72, cex.axis = 0.72, cex.lab = 0.72,
    mar = c(4.2,4,1,1) + .1)
plot(x = NULL, xlim = c(-4, 4), ylim = c(-3.5, 3), axes = FALSE,
     xlab = "Amount of confounding ({\\it α})", main = "Democrat Non-incumbents",
     ylab = "Estimated effect")
axis(side = 1, ps = 5)
axis(side = 2, las = 2)
abline(v = 0, col = rgb(0.1,0.1,0.1, alpha = 0.25))
polygon(x = c(confound[noninc.pos],rev(confound[noninc.pos])), y = c(noninc[noninc.pos, 3], rev(noninc[noninc.pos, 4])) , col = rgb(0.05, 0.05, 0.05,
                                                                                                                            alpha = 0.25), border =
        FALSE)
polygon(x = c(confound[noninc$non],rev(confound[noninc$non])), y = c(noninc[noninc$non, 3], rev(noninc[noninc$non, 4])) , col = rgb(0.05, 0.05, 0.05,
                                                                                                                            alpha = 0.25), border =
        FALSE, density = 25)
polygon(x = c(confound[noninc.neg],rev(confound[noninc.neg])), y = c(noninc[noninc.neg, 3], rev(noninc[noninc.neg, 4])) , col = rgb(0.5, 0.5, 0.5,
                                                                                                                            alpha = 0.25), border =
        FALSE)
lines(x = confound, y = noninc[,1], col = rgb(0.1,0.1,0.1),
      lwd = 2, lty = 2)
abline(h = 0, col = rgb(0.1, 0.1, 0.1, alpha = .25), lwd = 1)

text(x = 0, y = -3.5, "Negative campaigns\n worse off", pos = 2,
     cex = 0.5, col = rgb(0.5, 0.5, 0.5))

text(x = 0, y = -3.5, "Negative campaigns\n better off", pos = 4,
     cex = 0.5, col = rgb(0.5, 0.5, 0.5))
arrows(x0 = c(-0.2, 0.2), x1 = c(-2, 2), y0 = c(-2.7, -2.7), y1 =
       c(-2.7,-2.7), col = rgb(0.1, 0.1, 0.1, alpha = .25), length = 0.04)
plot(x = NULL, xlim = c(-4, 4), ylim = c(-3.5, 3), axes = FALSE,
     xlab = "Amount of confounding ({\\it α})", main = "Democract Incumbents",
     ylab = "Estimated effect", cex = .72)
axis(side = 1)
axis(side = 2, las = 2)
abline(v = 0, col = rgb(0.1,0.1,0.1, alpha = 0.25))
polygon(x = c(confound[inc.pos],rev(confound[inc.pos])), y = c(inc[inc.pos, 3], rev(inc[inc.pos, 4])) , col = rgb(0.05, 0.05, 0.05,
                                                                                                          alpha = 0.25), border =
        FALSE)
polygon(x = c(confound[1:22],rev(confound[1:22])), y = c(inc[1:22, 3], rev(inc[1:22, 4])) , col = rgb(0.05, 0.05, 0.05,
                                                                                              alpha = 0.25), border =
        FALSE, density = 25)
## polygon(x = c(confound[29:41],rev(confound[29:41])), y = c(inc[29:41, 3], rev(inc[29:41, 4])) , col = rgb(0.05, 0.05, 0.05,
##                                                                                                   alpha = 0.25), border =
##         FALSE, density = 25)

polygon(x = c(confound[inc.neg],rev(confound[inc.neg])), y = c(inc[inc.neg, 3], rev(inc[inc.neg, 4])) , col = rgb(0.5, 0.5, 0.5,
                                                                                                          alpha = 0.25), border =
        FALSE)
lines(x = confound, y = inc[,1], col = rgb(0.1,0.1,0.1),
      lwd = 2, lty = 2)
abline(h = 0, col = rgb(0.1, 0.1, 0.1, alpha = .25), lwd = 1)

text(x = 0, y = -3.5, "Negative campaigns\n worse off", pos = 2,
     cex = 0.5, col = rgb(0.5, 0.5, 0.5))

text(x = 0, y = -3.5, "Negative campaigns\n better off", pos = 4,
     cex = 0.5, col = rgb(0.5, 0.5, 0.5))
arrows(x0 = c(-0.2, 0.2), x1 = c(-2, 2), y0 = c(-2.7, -2.7), y1 =
       c(-2.7,-2.7), col = rgb(0.1, 0.1, 0.1, alpha = .25), length =
       0.04)
